/*
 * To change this template, choose Tools | Templates
 * and open the template in the editor.
 */
package TweetRay.Backgrounds;

import Raytracing.Constants;
import processing.core.PVector;

/**
 *
 * @author Matt
 */
public class ClearSkyBackground extends Background {

    private float turbidity = 3;
    public PVector toSun;
    /**
     * Angle in radians from the zenith to the sun vector.
     */
    public float thetaS;
    /**
     * Angle in radians to the sun vector from the south pole.
     */
    public float phiS;

    /**
     * Return the color of the sky for an angle in the sky.
     * @param dir Sky direction to sample
     * @return int color of the sky. Will be tonemapped using an exponential operator
     */
    public ClearSkyBackground() {
        toSun = new PVector(0, 0, 1);
        setSunVector(.5f, 0);
    }

    public void getColor(PVector dir, PVector ret) {


        if (dir.y <= 0 || toSun.y < 1.f) {
            ret.set(0, 0, 0);
        }
            toSun.normalize();

            float theta2, theta3, T, T2, chi, zenith_Y, zenith_x, zenith_y;
            float[] perez_Y = new float[5];
            float[] perez_x = new float[5];
            float[] perez_y = new float[5];

            theta2 = thetaS * thetaS;
            theta3 = theta2 * thetaS;
            T = turbidity;
            T2 = turbidity * turbidity;

            chi = (4.f / 9.f - T / 120.f) * (Constants.PI - 2.f * thetaS);

            zenith_Y = (4.0453f * T - 4.9710f) * (float) Math.tan(chi) - .2155f * T + 2.4192f;
            zenith_Y *= 1000.f;  // conversion from kcd/m^2 to cd/m^2

            zenith_x =
                    (0.00165f * theta3 - 0.00375f * theta2 + 0.00209f * thetaS + 0.0f) * T2
                    + (-0.02903f * theta3 + 0.06377f * theta2 - 0.03202f * thetaS + 0.00394f) * T
                    + (0.11693f * theta3 - 0.21196f * theta2 + 0.06052f * thetaS + 0.25886f);

            zenith_y =
                    (0.00275f * theta3 - 0.00610f * theta2 + 0.00317f * thetaS + 0.0f) * T2
                    + (-0.04214f * theta3 + 0.08970f * theta2 - 0.04153f * thetaS + 0.00516f) * T
                    + (0.15346f * theta3 - 0.26756f * theta2 + 0.06670f * thetaS + 0.26688f);

            perez_Y[0] = 0.17872f * T - 1.46303f;
            perez_Y[1] = -0.35540f * T + 0.42749f;
            perez_Y[2] = -0.02266f * T + 5.32505f;
            perez_Y[3] = 0.12064f * T - 2.57705f;
            perez_Y[4] = -0.06696f * T + 0.37027f;

            perez_x[0] = -0.01925f * T - 0.25922f;
            perez_x[1] = -0.06651f * T + 0.00081f;
            perez_x[2] = -0.00041f * T + 0.21247f;
            perez_x[3] = -0.06409f * T - 0.89887f;
            perez_x[4] = -0.00325f * T + 0.04517f;

            perez_y[0] = -0.01669f * T - 0.26078f;
            perez_y[1] = -0.09495f * T + 0.00921f;
            perez_y[2] = -0.00792f * T + 0.21023f;
            perez_y[3] = -0.04405f * T - 1.65369f;
            perez_y[4] = -0.01092f * T + 0.05291f;


            float theta = (float) Math.acos(dir.y);
            float phi = (float) Math.atan2(dir.z, dir.x);




            float gamma = AngleBetween(theta, phi, thetaS, phiS);

            // Compute xyY values
            float x = perezFunction(perez_x, gamma, theta, thetaS, zenith_x);
            float y = perezFunction(perez_y, gamma, theta, thetaS, zenith_y);
            float Y = perezFunction(perez_Y, gamma, theta, thetaS, zenith_Y);

            float X = (x / y) * Y;
            float Z = ((1.f - x - y) / y) * Y;


            float[] rgb = new float[3];
            rgb[0] = 1.f - (float) Math.exp(-(1.f / 15000.f) * (3.240479f * X - 1.537150f * Y - 0.498535f * Z));
            rgb[1] = 1.f - (float) Math.exp(-(1.f / 15000.f) * (-0.969256f * X + 1.875992f * Y + 0.041556f * Z));
            rgb[2] = 1.f - (float) Math.exp(-(1.f / 15000.f) * (0.055648f * X - 0.204043f * Y + 1.057311f * Z));


            if (rgb[0] < 0) rgb[0] = 0;
            
            ret.set(rgb[0], rgb[0], rgb[0]);


        
    }
    

    private float AngleBetween(float thetav, float phiv, float theta, float phi) {
        float cospsi = (float) (Math.sin(thetav) * Math.sin(theta) * Math.cos(phi - phiv) + Math.cos(thetav) * Math.cos(theta));
        if (cospsi > 1) {
            return 0;
        }
        if (cospsi < -1) {
            return Constants.PI;
        }
        return (float) Math.acos(cospsi);
    }

    private float perezFunction(float[] coeffs, float gamma, float theta, float thetaS, float lvz) {
        float A = coeffs[0], B = coeffs[1], C = coeffs[2], D = coeffs[3], E = coeffs[4];

        float den = (float) ((1.f + A * Math.exp(B)) * (1.f + C * Math.exp(D * thetaS) + E * Math.cos(thetaS) * Math.cos(thetaS)));
        float num = (float) ((1.f + A * Math.exp(B / Math.cos(theta))) * (1.f + C * Math.exp(D * gamma) + E * Math.cos(gamma) * Math.cos(gamma)));
        return (lvz * num / den);

    }

    /**
     * Set the sun vector using explicit spherical angles
     * @param theta Angle from the zenith to the sun in radians. Must be 0-PI.
     * @param phi Angle between the sun vector and the south vector in radians.
     */
    public void setSunVector(float theta, float phi) {
        thetaS = theta;
        phiS = phi;

        toSun.x = (float) (Math.cos(phi) * Math.sin(theta));
        toSun.z = (float) (Math.sin(phi) * Math.sin(theta));
        toSun.y = (float) Math.cos(theta);

    }
    /*
     *From IES Lighting Handbook pg 361
     * also, see http://www.gamedev.net/community/forums/topic.asp?topic_id=100346
     */

    /**
     * Sets the position of the sun using an earth position and time.
     * @param latitude 
     * @param longitude 
     * @param julianDay 
     * @param timeOfDay 
     * @param standardMeridian 
     */
    public void setSunByTime(float latitude, float longitude, int julianDay, float timeOfDay, float standardMeridian) {
        float solarTime, solarDeclination, solarAltitude,
                opp, adj,
                solarAzimuth;
        //

        standardMeridian *= 15.f;


        solarTime = timeOfDay
                + (float) (0.170f * Math.sin(4.f * Constants.PI * (julianDay - 80.f) / 373.f) - 0.129f * Math.sin(2.f * Constants.PI * (julianDay - 8.f) / 355.f))
                + (standardMeridian - longitude) / 15.0f;
        solarDeclination = (float) (0.4093f * Math.sin(2.f * Constants.PI * (julianDay - 81.f) / 368.f));
        solarAltitude = (float) (Math.asin(Math.sin(radians(latitude)) * Math.sin(solarDeclination))
                - Math.cos(radians(latitude)) * Math.cos(solarDeclination) * Math.cos(Constants.PI * solarTime / 12.f));

        opp = -(float) (Math.cos(solarDeclination) * Math.sin(Constants.PI * solarTime / 12.f));
        adj = -(float) ((Math.cos(radians(latitude)) * Math.sin(solarDeclination))
                + Math.sin(radians(latitude)) * Math.cos(solarDeclination) * Math.cos(Constants.PI * solarTime / 12.f));
        solarAzimuth = (float) Math.atan2(opp, adj);

        phiS = -solarAzimuth;
        thetaS = Constants.PI / 2.0f - solarAltitude;

        setSunVector(thetaS, phiS);
    }

    private float radians(float x) {
        return x / 180.f * Constants.PI;

    }
}
